// -----------------------------------------------------------------------------
// M A R I T I M E   R E S E A R C H   I N S T I T U T E   N E T H E R L A N D S
// -----------------------------------------------------------------------------
// Program  : lin_wacu
// File     : lin_wacu.cpp
// -----------------------------------------------------------------------------
// Revisions:
// Februari 13, 2009  asicherer
// -----------------------------------------------------------------------------
// testbed for the interactive waves model
// -----------------------------------------------------------------------------

// note: as this is testbed-only code: the architecture is not to be taken as an example

#define REMAIN // leave display on after simulation

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <deque>
#include <sys/time.h> // gettimeofday
#include "mermaid.h"
#include "FloatUtils.h"
#include "Array2D.h"
#include "InteractiveWaves.h"
#include "ShipForces.h"
#include "FileReader.h"
#include "defines.h"
#include "PulseControl.h"
#include "Timer.h"
#include "SolverFactory.h"
#ifdef WITH_PLOT
#include "OsgPlotter.h"
#include "PlotterFactory.h"
#include "ControlWindow.h"
#include "ShipGeometry.h"
#endif
#ifdef TEST
//#include <google/profiler.h>
#endif

// colors, shapes and sizes for waypoint/path display
#define PATHPOINT_SIZE    7        // size of path spheres
#define PATHPOINT_COLOR  "orange"
#define PATHPOINT_SHAPE  "sphere"
#define STARTPOINT_COLOR "green"
#define STARTPOINT_SHAPE "capsule"

REAL Fext = 0;

static int steps; ///< number of simulation steps
static int skips; ///< visualizion per 'skips' steps
static REAL currenttime = 0.0; ///< time stamp
static REAL timestep; ///< duration of one step in seconds
static REAL timestddev; ///< standard deviation of time step variation
static REAL pstart; ///< time to start pressure pulse
static int Nx1, Nx2; ///< size in points in x1,x2 direction
static REAL dx1, dx2; ///< grid step size in x1,x2 direction
static REAL grid; ///< minimum of grid step sizes
static REAL g; ///< gravity
static bool running = true; // simulation running?

static InteractiveWaves *pInteractiveWaves = NULL; ///< pointer to waves computing engine
static PulseControl *pPulseControl = NULL; ///< pointer to pressure pulse control
#ifdef WITH_PLOT
static ControlWindow *pControlWindow = NULL; ///< pointer to control window
static Plotter *pPlotterElevation = NULL; ///< pointer to plotter window
static Plotter *pPlotterNetElevation = NULL; ///< pointer to plotter window
static Plotter *pPlotterIntElevation = NULL; ///< pointer to plotter window
static Plotter *pPlotterBathymetry = NULL; ///< pointer to plotter window
static bool shipVisible = true; ///< is msh ship visible or not
#endif

#ifdef WITH_PLOT

/// draws path information on top of waves visualization
static void drawPaths()
{
    OsgPlotter *pl = dynamic_cast<OsgPlotter*>(pPlotterElevation);
    if (!pl)
    {
        return;
    }
    int nrp = pPulseControl->getNumberOfPulses();
    for (int ship = 0; ship < nrp; ++ship)
    {
        // draw point(s) where pulse(s) will come from
        PulseControl::Point start, end;
        bool ok = pPulseControl->getStartVector(ship, &start, &end);
        if (ok)
        {
            pl->addObject("1", STARTPOINT_SHAPE, STARTPOINT_COLOR, PATHPOINT_SIZE,
                          start.x1, start.x2, 1,
                          end.x1, end.x2, 1);
            // draw paths
            const PulseControl::Waypoints &path = pPulseControl->getPath(ship);
            int nn = path.size();
            for (int i = 0; i < nn; ++i)
            {
                pl->addObject("1", PATHPOINT_SHAPE, PATHPOINT_COLOR, PATHPOINT_SIZE,
                              path[i].x1, path[i].x2, 1);
            }
        }
    }
}

// callbacks for ControlWindow

/// sets exaggerated z scale
void setZscale(float scale)
{
    OsgPlotter *pl = dynamic_cast<OsgPlotter*>(pPlotterElevation);
    if (pl) pl->setZscale(scale);
    pl = dynamic_cast<OsgPlotter*>(pPlotterNetElevation);
    if (pl) pl->setZscale(scale);
    pl = dynamic_cast<OsgPlotter*>(pPlotterIntElevation);
    if (pl)  pl->setZscale(scale);
    pl = dynamic_cast<OsgPlotter*>(pPlotterBathymetry);
    if (pl)  pl->setZscale(scale);
}

/// saves waves file for later initialisation
void saveWaves()
{
    pInteractiveWaves->save();
}

/// returns number of waypoints of specified pulse (callback of ControlWindow)
int getNumberOfWaypoints(int ship)
{
    int nrp = pPulseControl->getNumberOfPulses();
    if (ship >= 0 && ship <= nrp)
    {
        const PulseControl::Waypoints &waypoints = pPulseControl->getWaypoints(ship);
        return waypoints.size();
    }
    return 0;
}

/// pauses and continues simulation
void toggleSimulation()
{
    running = !running;
}

/// pauses and continues simulation
void toggleMovement()
{
    pPulseControl->toggleMovement();
}

/// lowers ships
void startLowering()
{
    pPulseControl->startLowering(currenttime);
}

/// switches path display on or off (callback of ControlWindow)
void togglePathVisible()
{
    OsgPlotter *pl = dynamic_cast<OsgPlotter*>(pPlotterElevation);
    if (pl)
    {
        if (pControlWindow->isPathVisible() && pl)
        {
            while (pl->removeObject("1"))
            {
            }

        }
        else
        {
            drawPaths();
        }
    }
}

/// switches ship display on or off (callback of ControlWindow)
void toggleShipVisible()
{
    shipVisible = !shipVisible;
}



#endif // WITH_PLOT

/// reads test file
static bool setup(const char *testfile, const char *solver, Array2D<REAL> **ph0, Array2D<SHORT> **pwet, REAL *ph0_min, REAL *ph0_max)
{
    std::string areafile;
    try
    {
        FileReader fr(testfile);
        do
        {
            std::string section = fr.ReadString();
            if (section == "[Test]")
            {
                steps = fr.ReadInt();
                skips = fr.ReadInt();
                areafile = fr.ReadString();
            }
            else if (section == "[TimeSteps]")
            {
                timestep = fr.ReadFloat();
                timestddev = fr.ReadFloat();
                pstart = fr.ReadFloat();
                pInteractiveWaves = new InteractiveWaves();
                if (!pInteractiveWaves->setup(areafile.c_str(), solver, timestep, pstart))
                {
                    fprintf(stderr, "error setup %s\n", areafile.c_str());
                }
                pInteractiveWaves->getArea(&Nx1, &Nx2, &dx1, &dx2, &g);
                grid = std::min(dx1, dx2);
                pInteractiveWaves->getBathymetry(ph0, ph0_min, ph0_max, pwet);

#ifdef TEST
                printf("Matrix size: %d x %d\n", Nx1, Nx2);
#endif
            }
            else if (section == "[PressurePulse]")
            {
                ASSERT(Nx1 && Nx2, "");
                ASSERT(!pPulseControl, "");
                pPulseControl = new PulseControl(steps, timestep, Nx1, Nx2, dx1, dx2, pstart);
                if (!pPulseControl->setupPulses(fr))
                {
                    fprintf(stderr, "pulses not set up correctly\n");
                }
            }
            else if (section == "[PulsePaths]")
            {
                ASSERT(pPulseControl, "");
                pPulseControl->setupPaths(fr);
            }
        }
        while (!fr.Eof());
        return true;
    }
    catch (...)
    {
        return false;
    }
}


/// test program
int main(int argc, char *argv[])
{
    if (argc < 4)
    {
        std::string solvers = SolverFactory::getSolvers();
        fprintf(stderr, "lin_wacu testfileName solverName pg|osg|none [capture(0|1)] [realtime (0|1)]\n");
        fprintf(stderr, "solverName may be one of the following: ");
        fprintf(stderr, "%s\n", solvers.c_str());
        exit(-1);
    }
    const char *testfile = argv[1];
    const char *solver = argv[2];
#ifdef WITH_PLOT
    const char *plotter = argv[3];
#endif
    bool realtime = false;
    bool capture = false;
    capture = argc >= 5 && !strcmp(argv[4], "1");
    realtime = argc == 6 && !strcmp(argv[5], "1");
    if (capture)
    {
        puts("capture");
    }
    if (realtime)
    {
        puts("realtime");
    }

#ifdef USE_CUDAWAVESCOMPUTER
    puts("with CudaWavesComputer");
#else
    puts("normal WavesComputer");
#endif

#ifdef USE_CUDAINCOMINGWAVES
    puts("with CudaIncomingWaves");
#else
    puts("normal IncomingWaves");
#endif
    // setup test and bathymetry data
    Array2D<float> *ph0;
    Array2D<short> *pwet;
    REAL h0_min, h0_max;

    if (!setup(testfile, solver, &ph0, &pwet, &h0_min, &h0_max))
    {

        fprintf(stderr, "error setup %s\n", testfile);
    }

    // setup Froude-Krylov module (for now only for first ship)
    mml::Vector3d totalForce;
    mml::Vector3d totalTorque;
    ShipForces *pForces = NULL;
    const ShipGeometry* pShip = NULL;
    if (pPulseControl->getNumberOfPulses())
    {
        pShip = dynamic_cast<const ShipGeometry*>(pPulseControl->getGeometry(0));
        if (pShip)
        {
            pForces = new ShipForces (*pShip);
        }
    }

    // array for elevation data and pressure pulse
    Array2D<REAL> *pPressurePulse;
    Array2D<REAL> *pHeight = 0;
    Array2D<REAL> *pNetHeight = 0;
    Array2D<REAL> *pSpeed = 0;
    int nrp = pPulseControl->getNumberOfPulses();

#ifdef WITH_PLOT
 std::cout << "ENtered" ;
    // arrays for elevation data
    const int columns = Nx1 + 3;
    const int rows = Nx2 + 3;
    Array2D<float> waveheight(columns, rows);
    Array2D<float> netwaveheight(columns, rows);
    Array2D<float> bathymetry(columns, rows);
    Array2D<float> wetanddry(columns, rows);
    std::string bathymetryCaption(std::string(testfile) + " bathymetry");
    std::string wetanddryCaption(std::string(testfile) + " wet and dry points");
    std::string waveheightCaption(std::string(testfile) + " elevation");
    std::string netwaveheightCaption(std::string(testfile) + " corrected elevation");

    // create control window
    pControlWindow = new ControlWindow(&argc, &argv, nrp);

    // create bathymetry plotter window
    REAL pm = 10;
    REAL FMIN = -1.2 * fabs(pm);
    REAL FMAX = 1.2 * fabs(pm);
    if (h0_min == h0_max)
    {
        for (unsigned int i = 0; i < ph0->dim1(); ++i)
        {
            for (unsigned int j = 0; j < ph0->dim2(); ++j)
            {
                bathymetry[i][j] = (*ph0)[i][j];
            }
        }
    }
    else
    {
        for (int i2 = 1; i2 <= Nx2; ++i2)
        {
            for (int i1 = 1; i1 <= Nx1; ++i1)
            {
                bathymetry[i1][i2] = (h0_max - (*ph0)[i1][i2]) / (h0_max - h0_min) * (FMAX - FMIN) + FMIN;
            }
        }
    }
    pPlotterBathymetry = PlotterFactory::create(plotter, 1, 1, Nx1, Nx2, dx1, dx2, bathymetry, bathymetryCaption, FMIN, FMAX);

    // create wet and dry points plotter window
    for (int i2 = 1; i2 <= Nx2; ++i2)
    {
        for (int i1 = 1; i1 <= Nx1; ++i1)
        {
            wetanddry[i1][i2] = -10 * (*pwet)[i1][i2] + h0_min;
        }
    }
    Plotter *pPlotterWetDry = PlotterFactory::create(plotter, 1, 1, Nx1, Nx2, dx1, dx2, wetanddry, wetanddryCaption, FMIN, FMAX);

    // create net wave elevation plotter window
    netwaveheight = 0.0;
    pPlotterNetElevation = PlotterFactory::create(plotter, 1, 1, Nx1, Nx2, dx1, dx2, netwaveheight, netwaveheightCaption, FMIN, FMAX);
    OsgPlotter *pl = dynamic_cast<OsgPlotter*>(pPlotterNetElevation);
    if (pl)
    {
        pl->setZscale(1);
    }

    // create wave elevation plotter window
    waveheight = 0.0;
    pPlotterElevation = PlotterFactory::create(plotter, 1, 1, Nx1, Nx2, dx1, dx2, waveheight, waveheightCaption, FMIN, FMAX);
    pl = dynamic_cast<OsgPlotter*>(pPlotterElevation);
    if (pl)
    {
        pl->setZscale(1);
        pl->addAxes(osg::Vec3(1, 1, 1), 100, 2);
        pl->setDistance(1.7 * std::max(Nx1*dx1, Nx2*dx2));
    }

#endif

    // lower pulses automatically if there are only ShipPulses inside grid or MathPulses anywhere
    bool autostart = true;
    for (int i = 0; i < nrp; ++i)
    {
        const ShipGeometry* pShip = dynamic_cast<const ShipGeometry*>(pPulseControl->getGeometry(i));
        if (pShip) // msh ship
        {
            const mml::Vector3d& pos = pShip->getPosition();
            if (pos(0) < 0 || pos(0) > Nx1*dx1 || pos(1) < 0 || pos(1) > Nx2*dx2) // outside grid
            {
                autostart = false;
                break;
            }
        }
    }

    // start timing
    if (realtime)
    {
        getchar();
    }
    int m = 1;
    Timer timer;
#ifdef TEST
    double totalTime = 0.0;
    double totalWaveTime = 0.0;
    double totalCgTime = 0.0;
    double totalPulseTime = 0.0;
    double totalMoveTime = 0.0;
#endif

#ifdef TEST
//    ProfilerStart("lin_wacu.prof");
#endif

            if (autostart)

            {
                pPulseControl->startLowering(0.00000000001);
                autostart = false;
            }

    // Compute the wave evolution
    while (!steps || m <= steps)
    {
#ifdef WITH_PLOT
        pControlWindow->handleEvents();
#endif

        if (running)
        {

            double cgTime = 0.0;

#ifdef TEST
            double stepTime = seconds();
#endif
            // compute current time
            REAL dt = timestep + timestddev * FloatUtils::randomNumber(true); // new time step .. // AKASH: timestdev = 0 (No deviation, dt = constant)
            //    REAL dt = timestep * (FloatUtils::randomNumber(false) < 0.995 ? 1 : 2);
            currenttime = pInteractiveWaves->getTime() + dt;

            // compute pressure pulse
            double moveTime = seconds(); // AKASH: Gives the current time elapsed.
            pPulseControl->movePulses(currenttime);
            moveTime = seconds() - moveTime; // AKASH: This computes the time required for the moving a pulse.
#ifdef TEST
            totalMoveTime += moveTime;
#endif
            double pulseTime =seconds();
            pPulseControl->computePulsePattern(currenttime, 0.0, NULL, &pPressurePulse);
            pulseTime = seconds() - pulseTime;// AKASH: This computes the time required for computing Pulse Pattern
#ifdef TEST
            totalPulseTime += pulseTime;
#endif

            // compute new wave field
            double waveTime =seconds();
            pInteractiveWaves->compute(dt, *pPressurePulse, &pHeight, &pNetHeight, &pSpeed, &cgTime); // AKASH: What is pSpeed ? Ans: resulting speed potential
           // puts("after one compute");
            waveTime = seconds() - waveTime;
//            printf("\r waves:%.1f pulse:%.1f ms ", waveTime * 1000, pulseTime * 1000);
#ifdef TEST
            totalWaveTime += waveTime;
            totalTime += seconds() - stepTime;
            totalCgTime += cgTime;
#endif

            //        if (currenttime >= 99.98) exit(0);

            // compute Froude-Krylov force if the first pulse is a panel ship
            if (pForces)
            {
                pForces->computeFroudeKrylov(*pNetHeight, dx1, dx2, &totalForce, &totalTorque);
            }

            // plot the wave surface
            if (m % skips == 0)
            {
                /*REAL */Fext = 0.0;
                for (int i2 = 1; i2 <= Nx2; ++i2)
                {
                    for (int i1 = 1; i1 <= Nx1; ++i1)
                    {
                        Fext = std::max(Fext, (REAL)fabs((*pHeight)[i1][i2]));
    #ifdef WITH_PLOT
                        // tijdelijk tot hu file vervalt:
                        // show dry land as +4 meter
                        REAL depth = -(*ph0)[i1][i2];
                        if (!depth)
                        {
                            if (!(*ph0)[i1 - 1][i2] || !(*ph0)[i1 + 1][i2] ||
                                !(*ph0)[i1][i2 - 1] || !(*ph0)[i1][i2 + 1])
                            {
                                depth = 4;
                            }
                        } // !depth
                        waveheight[i1][i2] = std::max((*pHeight)[i1][i2], depth);
                        netwaveheight[i1][i2] = (*pNetHeight)[i1][i2];
    #endif
                    } // i1
                } // i2
    #ifdef TEST
               // printf("Fext = %10.5f\n", Fext);
    #endif
    #ifdef WITH_PLOT
                if (realtime)
                {
                    // slow down to real time
                    bool fast = false;
                    while (timer.getElapsedTime() < currenttime)
                    {
                        fast = true;
                    }
                    if (!fast)
                    {
                        puts("too slow");
                    }
                } // realtime
                char str[1000];
                if (pForces)
                {
                    sprintf(str, "h=%.1f dt=%.2f dx=%.2f step=%4d t=%4.2fs   fx=%.0f fy=%.0f fz=%.0f kN  mx=%.0f my=%.0f mz=%.0f kNm\n",
                        h0_max, timestep, grid, m, currenttime,
                        totalForce(0)/1000, totalForce(1)/1000, totalForce(2)/1000,
                        totalTorque(0)/1000, totalTorque(1)/1000, totalTorque(2)/1000);
                }
                else
                {
                    sprintf(str, "h=%.1f dt=%.2f dx=%.2f step=%4d; t=%4.2f\n",
                                h0_max, timestep, grid, m, currenttime);
                } // !pForces
//                printf("%s\n",str);
                waveheightCaption = str;

                // draw ship(s)
                for (int i = 0; i < nrp; ++i)
                {
                    OsgPlotter *pl = dynamic_cast<OsgPlotter*>(pPlotterElevation);
                    if (pl) // OpenSceneGraph plotter
                    {
                        // add ship to plotter
                        const ShipGeometry* pShip = dynamic_cast<const ShipGeometry*>(pPulseControl->getGeometry(i));
                        if (pShip)
                        {
                            char buffer[10];
                            sprintf(buffer, "ship%02d", i);
                            pl->removeObject(buffer);
                            if (shipVisible)
                            {
                                std::vector<ShipGeometry::ExtendedPanel> panels;
                                pShip->getExtendedPanels(&panels);
                                int nPanels = pShip->getNumberOfPanels();
                                pl->addPanelRoot(buffer);
                                for (int i = 0; i < nPanels; ++i)
                                {
                                    ShipGeometry::ExtendedPanel& panel = panels[i];
                                    pl->addPanel(buffer, "yellow",
                                                 panel.vertices, panel.normal);
                                } // for all panels
                            } // shipVisible
                        } // pShip
                    } // pl
                    if (pControlWindow->isPathVisible() && !pl /* no OpenSceneGraph plotter */)
                    {
                        const PulseControl::Waypoints &path = pPulseControl->getPath(i);
                        unsigned int nn = path.size();
                        for (unsigned int w = 0; w < nn; ++w)
                        {
                            waveheight[int(path[w].x1 / dx1 + 0.5)][int(path[w].x2 / dx2 + 0.5)] = 10;
                        }
                    } // pathVisible
                } // for all ships
                pPlotterElevation->plot();
    #endif // plot
            } // if (m % skips == 0)
            ++m;

#ifdef delft
            if (m>steps-9)
            {
                int ccc = 9-(steps-m);
                printf("%2d %15.12f\n", ccc, cgTime);
                char buffer[20];
                sprintf(buffer,"m_bpsi_rl_%02d.txt",ccc);
                rename("m_bpsi_rl.txt", buffer);
                sprintf(buffer,"m_psi_sol_%02d.txt",ccc);
                rename("m_psi_sol.txt", buffer);
            }
#endif
        } // running
    } // for all steps or indefinitely

#ifdef TEST
//    ProfilerStop();
    printf("Computation   : %G seconds per step\n", totalTime / steps);
    printf("Move pulses   : %G seconds per step (%4.1f%% of total computation)\n", totalMoveTime / steps, totalMoveTime / totalTime * 100);
    printf("ShipPulse     : %G seconds per step (%4.1f%% of total computation)\n", totalPulseTime / steps, totalPulseTime / totalTime * 100);
    printf("Waves         : %G seconds per step (%4.1f%% of total computation)\n", totalWaveTime / steps, totalWaveTime / totalTime * 100);
    printf("with Solver   : %G seconds per step (%4.1f%% of waves computation)\n", totalCgTime / steps, totalCgTime / totalWaveTime * 100);
    printf("Performance   : %G seconds per step per element\n", totalTime / steps / (Nx1 * Nx2));

#endif

#ifdef WITH_PLOT
    if (capture)
    {
        pl->toggleScreenCast();
        pl->plot();
//        pl->setZview();
  //      pl->plot();
    }
 #ifdef REMAIN
    for (;;)
    {
        pControlWindow->handleEvents();
        pPlotterElevation->plot();
        pPlotterNetElevation->plot();
    }
#endif
    delete pPlotterBathymetry;
    delete pPlotterWetDry;
    delete pPlotterElevation;
    delete pPlotterNetElevation;
    delete pControlWindow;
#endif
    delete pInteractiveWaves;
    delete pPulseControl;
    delete pForces;

    return 0;
}

